Renewable and Sustainable Energy Reviews 16 (2012) 1709-1720 


Renewable and Sustainable Energy Reviews 


journal homepage: www.elsevier.com/locate/rser 


Contents lists available at SciVerse ScienceDirect 


Characterization of solar flat plate collectors 


F. Cruz-Peragon, J.M. Palomar, P.J. Casanova”, M.P. Dorado‘, F. Manzano-Agugliaro4%* 


4 Dept. of Mechanics and Mining Engineering, EPS de Jaen, Universidad de Jaen, Campus Las Lagunillas s/n, 23071 Jaen, Spain 
> Dept. of Electronics Engineering and Automation, ESP de Jaen, Universidad de Jaen, Campus Las Lagunillas s/n, 23071 Jaen, Spain 
€ Dept. of Physical Chemistry and Applied Thermodynamics, EPS, Campus de Rabanales, Universidad de Cordoba, Campus de Excelencia Internacional Agroalimentario, 


ceiA3, 14071 Cordoba, Spain 


4 Dept. Rural Engineering, University of Almeria, 04120 Almeria, Spain 


ARTICLE INFO 


Article history: 

Received 8 March 2011 

Received in revised form 

21 November 2011 

Accepted 25 November 2011 
Available online 18 January 2012 


Keywords: 

Model optimization 
Numerical approach 
Solar collector 
Correlation functions 


Contents 
1. Introduction and previous work 
2. Materials and methods 


2.1. Equipment.............. 
2.2. Collector model 
2.2.1. Steady state equations 


ABSTRACT 


Characterization of solar collectors is based on experimental techniques next to validation of associ- 
ated models. Both techniques may be adopted assuming different complexities. In this work, a general 
methodology to validate a collector model, with undetermined associated complexity, is presented. It 
serves to characterize the device by means of critical coefficients, such as the film (convection) transfer 
coefficient, plate absortance or emmitance. The first step consists of identifying those significant param- 
eters that match the selected model with the experimental data, via nonlinear optimization techniques, 
applied to steady state conditions. Second, new correlations must be adopted, in those terms where it is 
necessary (i.e. film coefficient equations). Finally, the overall model must be checked in transient regime. 
To illustrate the technique, a tailor-made prototype flat plate solar collector has been analyzed. An inter- 
mediate complex collector model has been proposed (2D finite-difference method). Both steady and 
transient states were analyzed under different operating conditions. Parameter identification is based on 
Newton's method optimization. For parameter approximation, exponential regression functions through 
multivariate analysis of variance is proposed among many other alternatives. Results depicted a robust- 
ness of the overall proposed method as starting point to optimize models applied to solar collectors. 


© 2011 Elsevier Ltd. All rights reserved. 
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1. Introduction and previous work 


Processes of industrialization and economic development 
require important energy inputs [1]. Fuels are the world’s main 
energy resource and are considered the center of energy demands. 
However, reserves of fossil fuels are limited and their large-scale 
use is associated with environmental deterioration [2]. These facts 
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Nomenclature 
Tfp fluid temperature into the upper and lower horizon- 
a constant factor of the regression function T tal pipes(K) 
but exponent (constant value) of the regression function Tp material temperature of the upper and lower hori- 
T zontal pipes (K) 
Bi Biot number Tm mean top surface temperature (K) 
Cf fluid specific heat (J/(kg K)) Tw mean temperature in the inner side of the pipe in a 
Cp specific heat of the absorber plate (J/(kg K)) differential element (K) 
Di inner diameter of the risers (m) UP" Unknown parameters vector at a stationary state at 
e plate thickness (m) time ‘n’ 
eq 2D equivalent thickness of the pipe-weld-plate (m) UL volumetric heat transfer loss coefficient of the col- 
er pipe thickness (m) lector (W/(m? K)) 
fi exponential regression function ‘I’ U; global heat loss coefficient of the collector 
Gr Grashof number (W/(m? K)) 
h final estimated film coefficient of the inner wall Up heat loss behind the collector (W/(m? K)) 
water (W/(m? K)) U: top heat loss (W/(m? K)) 
ho initial film coefficient of the inner wall-water Uc volumetric equivalent heat transfer coefficient of 
(W/(m? K)) the control volume pipe-weld-plate (W/(m? K)) 
Ir global incident radiation over tilted surface (W/m?) Xul standardized input parameter ‘u’ of a regression 
ke extinction coefficient of the covering glass (m~!) function for output parameter ‘I’ 
kf pipe conductance (W/(m K)) x transversal direction distance in the collector (m) 
kp conductance of the absorber material (W/(m K)) y longitudinal direction distance in the collector (m) 
kh proportionality factor applied to h a thermal diffusivity of the absorber plate material 
kq proportionality parameter associated to the volu- (copper) (m?/s) 
metric equivalent heat transfer coefficient Uc ag directional absorbance of the cooper black plate 
kuc proportionality factor applied to U; A increment 
kir conduction heat transfer coefficient of the insula- Eg glass emmitance 
tion material (W/(m K)) Ep plate emmitance 
Ke global heat transmission coefficient plate-fluid p absorber density (kg/m?) 
(W/(m2 K)) Pf refrigerant density (kg/m?) 
my fluid flow rate into the risers (kg/s) o latitude (rad) 
i, fluid flow rate in lower and upper pipes (kg/s) w timing angle (rad) 
n2 refraction index of the covering glass ô declination (rad) 
Nu Nusselt number y azimutal angle (rad) 
Pr Prandtl number 0 incident angle of beam irradiance over the tilted sur- 
qi unitary surface heat radiation on the absorber face (rad) 
(W/m?) y objective function to minimize by the identification 
qı, S internal heat generation in the absorber due to solar process 
radiation (W/m?) r auxiliary expressions for Eqs. (7.a), (7.b), (9.a) and 
q'2 unitary surface heat loss in the absorber plate (9.b) 
(W/m?) 
q2 loss of the internal heat generation in the absorber Sub-index 
(W/m?) i x direction 
q'3 unitary surface heat transferred to the working fluid j y direction 
(W/m?) m number of nodes in x direction 
a3 loss of the internal heat generation in the absorber w number of input parameters in a regression function 
due to heat transmission to water (W/m?) P number of riser 
dv internal heat generation in the absorber (W/m?) r node number 
q'u net heat released in the refrigerant water per length 
unit (W/m) Super-index 
fi inner radius of the riser (m) k time in transient regime analysis 
Te outer radius of the riser (m) n interval of time 
T'i inner radius of lower and upper pipe (m) 
Re Reynolds number 
—_ collector tilt (rad) have encouraged growth in the use of renewable energy resources 
SP known conditions vector at time ‘n’ worldwide [3]. Solar energy is considered one of the main promis- 
t time (s) ing alternative sources of energy to replace the dependency on 
T temperature (K) other fossil energy resources [4,5]. Solar water heating systems 
Ta cuidado si es To (Eq. (2.b)) are very common systems, extensively used in many countries 
To room temperature (K) with high solar radiation potential, such as Mediterranean coun- 
Tp plate temperature in the bulb thermometer (K) tries [6]. They are often viable to replace fossil fuels used for many 
Ta refrigerant fluid inlet temperature (K) home applications [7]. Conventional flat plate solar collectors with 
Tyo refrigerant fluid outlet temperature (K) a metal absorber plate and covers are used to transform solar energy 
Ty fluid mean temperature in a differential element (K) into heat [8]. Optimization of solar heating systems provides both 


running and design parameters that ensure maximal collected heat. 


F. Cruz-Peragon et al. / Renewable and Sustainable Energy Reviews 16 (2012) 1709-1720 1711 


Main parameters to be considered are related with both geometry 
and materials properties [9-16]. An optimal design must ensure 
minimal working fluid pressure loss, besides a nearly constant 
temperature on each transversal section of the plate. Otherwise, 
temperature between pipes will rise, thus increasing radiation heat 
loss. To carry out the optimization of the process, a system sim- 
ulation, including model validation and analysis under different 
working conditions are needed. In this sense, the first step consists 
of determining the parameters associated to heat collection and 
transmission. Subsequently, a mathematical model considering the 
following premises must be established: 


a) Parameters associated to the physical properties of the system. 
According to the literature, main material properties that influ- 
ence the system are well known [17-19]. Conventional absorber 
plates usually exhibit either parallel or serial hydraulic con- 
figurations. These configurations (including “Z” configuration) 
besides circulating regimes showing small Reynolds number 
(Re) lead to non-uniform flow distribution. On the other hand, 
heat transfer from wall to working fluid leads to non-uniform 
transversal temperature distribution on the plate. Thus, disre- 
gard of some critical parameters related to both fluid flow and 
heat transfer by convection makes difficult the optimal design of 
the system. Although primary surface pressure loss can be easily 
estimated [20,21], little research about coefficients of pressure 
loss in junctions and bifurcations in collector pipes has been 
done. In fact, previous works only consider parameters depen- 
dent on Re at different flow ratios [22,23]. Flow ratio is defined 
as the relation between the discharge in each riser (being riser 
no. 1 the nearest to the inlet port) and the total discharge of the 
collector. In this sense, there is an extensive literature devoted 
to semi-empirical correlations for film coefficients approxima- 
tion [24-26]; however, correlations for pressure loss in junctions 
and bifurcations are still needed [27,28]. On the other hand, 
extremely small Re means axial flow velocity is distorted. This 
feature makes difficult to understand and approximate the pro- 
cess [29]. Also, flow distribution seems to be strongly influenced 
by the Grashof number (Gr) and the pipe leans [30,31]. 

b) Simulation complexity. Collector characteristics are usually well 
known and needed for model design. However, system char- 
acterization may increase the complexity of the model. This 
complexity can be described under both stationary and transient 
regimes through computational fluid dynamics (CFD) methods 
[19,32-35]. The target is to reach a compromise between the 
complexity of the selected model and the accuracy of its results, 
which should adequately simulate the system behavior. 


Single-dimensional methods are able to simulate with high accu- 
racy the majority of the situations. Nevertheless, to carry out an 
optimal design, two-dimensional methods based on finite differ- 
ence or control volume techniques must be considered, so a good 
fitting can be achieved [36,37]. In this sense, to determine the tem- 
perature distribution over the plate, an adequate simulation of 
the heating system entails an associated mathematical difficulty, 
considering both characterization and CFD techniques. In fact, in 
case experimental data show high temperature difference between 
upper tubes, previous models are useless [36,38,39]. 

For these reasons, validation of a defined numerical or math- 
ematical model involves data acquisition and identification of 
parameters that show abnormal behavior according to the liter- 
ature. Eventually, it can force to re-define some of the equations 
associated to the model. 

Thus, the objective of this work is to propose a complete pro- 
cedure to optimize the equations associated to a collector model. 
In the first step, the identification of the parameters associated to 
materials and heat transfer provides the system validation over 


the time. Then, previous values will help to provide new corre- 
lation functions associated to the collector behavior. In sum, a first 
approach to optimize a two-dimensional finite difference model 
associated to a flat plate solar collector with low flow rate is pro- 
posed. 

Identification of a thermo-physical process involves the deter- 
mination of its parameters by solving an inverse heat transfer 
problem (IHTP). Inverse problems consist in the determination 
of the initial parameters when solutions are known. This kind of 
analysis has been proposed by many researchers during the past 
years [40,41]. Also, it has been proved that solutions to IHTP are 
usually unique, though estimations are not always numerically 
stable [42]. In this work, an inverse problem (identification pro- 
cedure) is performed by means of iterations over a direct problem 
(system model) using a derivative dependent method. The identifi- 
cation is based on a non-linear optimization technique that uses an 
objective function, Y, obtained from the sum of squares of the dif- 
ferences between the model and the experimental values [43,44]. 
This quadratic form provides the best performance in a wide range 
of engineering optimization problems. This technique applies the 
second order Newton minimization method, in combination with 
a Taylor expansion series of W, truncated to the third derivative 
[43]. According to the literature, it is observed that the solution set 
(temperatures) is continuous and smooth, without drastic changes 
[36,38,39]. The same assumptions have been applied to heat trans- 
fer problems [45]. For these reasons, when a quasi-Newton method 
truncated to the second derivative is selected, the computational 
efforts are reduced while keeping their accuracy [43]. If the num- 
ber of measured points (in terms of time-space) is not lesser 
than the number of parameters to be identified, the solution of 
a non-linear steady IHTP is unique. The simultaneous identifica- 
tion of thermo-physical characteristics, internal heat source and 
heat flux would be unique if either before or during the numer- 
ical experiment, the monotonicity or piecewise-monotonicity of 
the parameters to be determined is found to be likewise [40]. In 
the present work, none of these premises fulfil the problem condi- 
tions because only two sensors were used. For this reason, results 
only describe a general approach. Nevertheless, specificity could be 
increased by adding more sensors over the absorber and into the 
risers. 


2. Materials and methods 
2.1. Equipment 


The solar collector system, as shown in Fig. 1, is composed ofa set 
of parallel copper pipes under a ‘Z’ disposition. They are attached to 
a copper flat plate (absorber) by means of a tin-silver (6%) welding. 
This assembly is covered with a matt black paint coating, operat- 
ing as a selective surface, thus making possible a great radiation 
collection. The casing is made of a white-lacquered aluminium 
profile. Two frames of this type are superposed. One of them con- 
tains white glass, while the other one contains the pipes-plate unit. 
The isolator is composed of two thin aluminium sheets separated 
by an expanded polystyrene layer. It incorporates inlet and outlet 
water thermometers and manometers, next to a bulb thermome- 
ter placed on the absorber plate. Security valves are attached to this 
equipment. The main characteristics and material properties of this 
collector are shown in Table 1. They provide an accurate simulation 
of the system [46]. 


2.2. Collector model 


Analyses must be carried out under rigorous model condi- 
tions, including the main features of a representative dynamic 
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Flow distribution 


LEGEND 

2.- Thermometer in inlet pipe 
4. Pipes and risers 

6.- White glass covering 

8. Purging valve 

10.- Thermometer in outlet pipe 
12.- Thermometer of flat plate 


1.- Manometer in inlet pipe 

3.- Cooper absorber plate 

5.- Isolators 

7- Aluminium case r 

9.- Manometer in oulet pipe 

11- Thermometer of glass covering 


Fig. 1. Collector under study. 


Table 1 

Collector properties. 
Property Value 
Absorber material Copper 
Absorber thickness (mm) 2 
Absorber useful surface (m?) 0.21 
Inner diameter of the horizontal tubes, R, (mm) 20 
Inner diameter of the riser, R; (mm) 10 
Number of risers 15 
Riser length (mm) 450 
Distance between risers (mm) 30 
Glazing thickness (mm) 4 
Complete useful surface (m?) 0.37 


Posterior insulation (mm) 10 (aluminium 
sandwich)+7 
(polyurethane) 


Collector plate efficiency factor, Fr 0.7 


simulation model and applied only to the plate-risers unit [47]. In 
this particular case, the following conditions were adopted: 


(i) A heat transfer analysis by finite difference approach of the 
absorber plate and the fluid that flows through the pipes must 
be carried out. Some simplifications consider the plate as a two- 
dimension system and each pipe as a one-dimension system. 
The model consists of a set of non linear equations for both 
steady and transient states, which must be solved by iterations 
[48]. 

(ii) Conduction coefficients and diffusivities of the plate, pipes 
and cover, as well as insulating materials must be taken into 
consideration. On the other hand, densities, specific heat, vis- 
cosity, conductivity and Prandtl number of the circulating fluid 
(water) are approximated by 4th-order polynomial functions, 


i Plate portion i 
4 > 


Nodesi: 12 34567 


12345 67 
Finite difference approach 


xt 
i Nodes i = 1 to 7 


Fig. 2. Sample network. 


depending on the fluid temperature [49]. The necessary data 
are collected from the literature [50]. 

(iii) Heat loss coefficients of the plate are considered to be depen- 
dent on both average and room temperatures, at each interval 
of time. 

(iv) Temporary variations of irradiation, room and water inlet tem- 
peratures are considered model inputs. In this sense, both the 
unitary heat release of the tilted surface (adjusted with the 
help of an isotropic model, with reflectance equal to 0.2) and 
the absorptance-transmittance product depend on time [17]. 
Additionally, time variations of water outlet temperature and 
bulb temperature (over the flat plate) are known but only used 
during the identification process. 

(v) According to Streeter et al., flow distribution into each pipe 
may be calculated by means of the friction loss and conven- 
tional techniques [20]. However, secondary loss coefficients 
applied to collector junctions and bifurcations are both Re (of 
the main pipe flow) and flow rate dependent [22]. According 
to Wietbrecht et al., they may be approximated by exponential 
functions [22]. 

(vi) Film heat transfer coefficients are initially approached by 
means of semi-empirical equations [24,25,28,50,51]. Although 
there are many particular correlations, those listed in Incropera 
et al. have been selected [50]. Once the optimization process is 
finished, the coefficients are empirically evaluated to validate 
the model. 


To carry out the analysis, the absorber surface was divided into 
multiple similar units, using a mesh, so it could be considered a dis- 
crete problem. A grid sample is shown in Fig. 2. According to this, 
the complete study could be approximated by means of finite differ- 
ence technique [52]. Because pipes distribution along the absorber 
plate is considered to be constant, only a portion of the plate was 
examined. The plate is attached to 15 identical tubes, symmetri- 
cally distributed regarding the y-axis, thus dividing the plate in 30 
similar zones. However, only half of the area between tubes is con- 
sidered, provided the symmetry in the geometric design. For this 
reason, only 14 similar zones were considered, as shown in Fig. 2, in 
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Fig. 3. Flow distribution ratio on each riser. 


addition to another two half zones in both plate extremes, regard- 
ing x-axis. The boundary line is considered adiabatic, due to the 
presence of the isolator. Boundary loss is around 3% of the net heat 
released [49], so it was neglected. 

Flow distribution along the risers can be considered constant 
unless the regime is laminar [22], which happens when the flow 
rate is very low [46,49]. In this work, Re is smaller than 100 along 
the circuit of water, due to low flow rate and non-developed flow, 
making viscous forces more important than inertia ones. Primary 
hydraulic surface loss in laminar flow can be calculated by the 
Poiseuille equation, but form loss (secondary loss) analysis needs 
the application of some experimental coefficients [20]. With the 
help of these coefficients, pressure loss and flow distribution along 
the pipes may be analyzed. Another secondary kind of coefficient 
takes extremely high values, for Re below 100. So, it makes pos- 
sible to calculate form loss by a general exponential expression 
depending on Re and the flow ratio in junctions. It derives from 
the empirical and validated functions observed by Weitbrecht et al. 
[22]. 

The Fourier law, related to heat conduction through a material, 
given by Eq. (1), describes the system dynamics as follows: 


w Ary, e TT PT PT w _ 10T 


PCp ðr ” ðx? ° dy2 Əz? kp wat 
; a 1 
with A = a (1) 


where T is the temperature (K), t is the instantaneous time (s), 
kp is the conduction coefficient of the absorber (W/(m K)), «œ is its 
thermal diffusivity (m/s), p is its density (kg/m?), cp is its specific 
heat value (J/(kgK)), and y is the internal heat generation in the 
absorber (W/m?). 

In fact, the dependent variables presented in Eq. (1) must be 
optimized by a method that will be explained in detail later. The 
flow rate differs along each riser, as shown in Fig. 3. The main 
consequence is the different temperature distribution between the 
sections of the plate. So, the developed model of heat transfer must 
be adapted to each riser attached to the flat plate. For this reason, 
the study is performed in 14 similar portions of pipes of the same 
geometry (including two half zones in both plate extremes), but 
with different water flow rate due to hydraulic unbalance caused 
by the high form loss coefficients associated to the laminar regime. 

The complete heat balance considers boundary conditions, that 
is, the input and output of energy to/from the plate at each time 
interval and per surface unit (W/m2). Firstly, the unitary surface 


heat radiation, qi (W/m2), is determined from the global incident 
radiation over the tilted absorber, Ir (W/m), reduced by the prod- 
uct transmittance—absorptance of the ‘ta’ glass-absorber unit, as 
shown in Eq. (2.a). 


qi(t) = Irta(t) (2.a) 


Secondly, the unitary heat loss of the plate, q% (W/m2), is calcu- 
lated from the temperature difference between the selected point 
and the environment, as follows (Eq. (2.b)): 


q3(t) = Ut (T(t) — Ta(t)) with U, = UAT(t), Ta(t)}, (2.b) 


The global heat loss coefficient, U; (W/(m? K)), is composed of two 
factors: (a) the loss behind the collector, Up (W/(m? K)), due to 
the heat conduction through the isolator and the casing, and (b) 
the convection loss to the atmosphere. Its values depend on the 
material properties and wind speed. The top loss, U; (W/(m? K)), 
originates from solar radiation and re-radiation, system tilt, orien- 
tation and plate temperature. 

Finally, q3” (W/m?) (associated to pipes geometry) provides the 
available heat that passes through the working fluid, named use- 
ful heat. It may be calculated by means of the film coefficient, h 
[W/(m? K)], related to the convection heat transfer from the inner 
wall of the pipe, with mean temperature Ty (K) and water mean 
temperature T; (K), along a differential length Ay (m) in direction 
y, as shown in Eq. (2.c). 


(pcp AT; + Aystr;? ppcp(dTs/dt)) 
(27r; Ay) 


q3(t) = h(t)(Tw(t) — T(t) = 


(2.c) 


where rhy is the fluid flow rate (kg/s), To is the room temperature (K), 
cris the water specific heat value (J/(kg K)), opis the density (kg/m), 
AT; is the increment of temperature (K) and r; is the internal radius 
of each riser (m). Finally, this procedure comprises some heat input 
to the water through the horizontal manifolds before and after the 
risers. For their determination, Eq. (3) is included. This equation 
provides the inlet temperature value into each riser. 

The very low water flow rate and the small scale provide a non- 
developed flow into the pipes. Re is small enough to distort the axial 
speed of the fluid along the riser by buoyancy forces [29]. The filed 
flow is strongly influenced by Gr and the elevation angle of the riser 
[30]. It makes impossible to assure the capability of the equations 
provided by the literature for the determination of the film coeffi- 
cient, h, only by CFD techniques. Provided the difficulty to carry out 
an initial approach to achieve an accurate value of h, it is calculated 
by means of an initial value, ho, obtained from semi-empirical rela- 
tions of heat transfer theory and multiplied by a proportionality 
term associated to the coefficient, kp, as shown in Eq. (3). To assure 
coherent results, the optimal value for the coefficient is needed. 


h(t) = 10**ho(t); ho = hofkp, ri, Nu} (3) 


On the other hand, when Biot number, Bi, is evaluated (Bi=he/kp, 
where e is the thickness of the plate), it gives a value much smaller 
than one. This means the temperature across the thickness of the 
plate can be considered constant, giving the same results in both 
sides of the plate, as found by other researchers [52]. So, the unit 
plate-pipes can be analyzed as a two-dimension system. Accord- 
ing to this, the internal heat generation of Eq. (1) results from the 
boundary conditions of Eq. (2), adapted to Eq. (4), as follows: 


a(t) = q(t) — gat) — 4a(t)in W/m?; qu(t)=S = ne. 
a= 2 ana u=” 1a) 


As can be seen, this simplification leads to the introduction of 
a volumetric heat transfer loss coefficient, U, (W/(m?K)), being 
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necessary to analyze in detail the useful heat. As a result of the 
proposed two-dimension problem, the temperature of the unit 
plate-pipes across its section is assumed to be constant. Concern- 
ing the risers, this assumption is not accurate because an inner wall 
temperature for each riser, Tw, is established. This temperature is 
different from the mean top surface temperature, Tm; temperature 
set varies from the top to the inner wall. The weld that connects 
the pipes to the plate is also included in the study. These consider- 
ations, besides the geometry and materials of each control volume, 
difficult the calculation of the global heat transfer coefficient of 
the plate-fluid, Kc [W/(m? K)]. So, the equation of useful heat col- 
lected by the working fluid must be recalculated. The simplification 
from three to two dimensions obliges to accommodate the prob- 
lem. For this purpose, an equivalent thickness, eg (m), in addition 
to a volumetric equivalent heat transfer coefficient, Uc [W/(m?K)], 
that considers the heat from the top to the inner wall, are evaluated. 
These terms must include a proportional coefficient depending on 
the pipe geometry, plate weld and position, called kyc. As happened 
to kp term, the optimal value that provides an adequate behavior 
of the collector depends on time. So the heat released to the fluid 
matches Eq. (5): 


q(t 
43(E) = KeltX Tm (t) = TC) a= 8, 
2_ r2 
e] = (e Ax + x(r — r?)). ee ee 


=j 
Te 1 Te Te . 7 
( AE)” kp n(; )) + 43(t) = U(t T(E) = TCE); 


-1 
Te 1 | fe Te 
(Ea ! en (£)) " 
Eqs. (4) and (5) help to determine the temperature in particular 
points of the plate. The one dimensional useful heat that passes 


to the water, q(t) (W/m), considering steady state conditions and 
Ay, may be calculated according to Eq. (6) as follows: 


TaD gy = Ry = Tw ~ Tp) (6) 


2.2.1. Steady state equations 

For each (i,j) point of the mesh, and depending on the increments 
Ax and Ay on the absorber fin, the numeric adaptation of the set 
of Eqs. (1)-(6) for steady state analysis follows the expressions: 


U; 


L 
general form : Ti-1 j + Tiga j + Tij-1 + Tij+1 — Tij 4 } Rep Ax as] — 


considering both limits (Toj and Tm+ı correspond to fluid tempera- 
tures associated to nodes r=1 and r=™m, respectively). 


Tijl t To j-1 Tint + Ting j—1 


27r; Ayh 
—ı;(1+T')=0; r= +; 
m 1,( ) Qing Cp 


(7.b) 
The set of (m+2)n Eqs. (7.a) and (7.b) is used to establish the basic 
model in the plate-pipes unit, considering the unknown tempera- 
tures on the plate (two-dimension) and the working fluid passing 
through the tube (one-dimension). This is referred to the system 
portion described above (Fig. 2). 

A non-linear set of equations that must be solved by iteration 
is presented. For this purpose, plate temperature is considered to 
be room temperature. Then, the rest of the parameters such as Uz, 
Uc, h and S will be calculated. A matrix calculus is built-up, result- 
ing in a new set of plate temperatures. This process continues until 
temperatures converge to constant values. It is very important to 
advise that this iterative process belongs to the optimization pro- 
cess explained below, considering known all the parameters that 
define UP". 

To evaluate the set of Eq. (7) for each riser, the inlet water tem- 
perature must be estimated. For this purpose, a light modification 
is adopted and applied to the input horizontal manifold. In this 
pipe, the temperature is evaluated for only one dimension (x-axis), 
between the junctions that connect two consecutive risers. Pro- 
vided that the global input flow temperature and flow distribution 
throughout the unit are known, the equations set can be success- 
fully solved. It gives information about the heat delivered to the 
water in the lower pipe (with inner radius 1’;), from the riser p to 
the riser p+ 1, separated Al. It results in the following expression: 


Ur + Uc , i Uc P j 
Tp-1 + Tp+1 — Tp |24 A Ax’ Ay tiop AX Ay 
Ax’ Ay’ 
aaa Y (S+ ULT) (8.a) 
p 
’ j . , 2ar, Alh. 
Tp pl” +Tfp-1 -Tip 0 +I")=0; I'= 2rn,Cp i (8.b) 


where Ax = Al, Ay'=2r'i. The fluid that flows through the lower 
pipe is non-uniform and continuously reduced, as shown in Fig. 
3. Once the evaluation of Eq. (7) is finished, the upper pipe is 
analyzed. This will provide the outlet global temperature of the 
fluid. 

The proposed numerical analysis has been compared to other 
analytical methods from the literature, providing similar results for 
the film coefficient h and same qualitative tendencies for pressure 
loss [53]. 


Ax Ay 


kp (S+ ULTo) 


adiabatic riser boundary form : 2Tj_4 j + Tij-1 + Tijy1 — Tij 4 Ur + Uc Ax ay! + Tia j - Ax Ay = 
p 


U, 


adiabatic fin boundary form : 2Ti41 j + Ti j-1 + Tij41 — Tij 4 H i 


For a finite different approach, m corresponds to the number of 
nodes in ‘x’ direction (i varies from 1 to m). Moreover, the work- 
ing fluid receives a certain quantity of heat when passing from y =j 
to j+1. It flows into the riser placed at each extreme of the mesh, 
as shown in Fig. 2. Thus, an additional equation for each one, i.e. 
Eq. (6), is needed. A numerical approach is presented in Eq. (7.b), 


AX as] = OV S4 ens 
P 


kp 


2.2.2. Transient state equations 

When Eqs. (7.a), (7.b), (8.a) and (8.b) are used under transient 
state condition, the inclusion of settings that consider the deriva- 
tive of temperature vs. time is required. To evaluate it, an implicit 
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method has been proposed. The method is very stable and does 
not require extremely narrow intervals of time. In this case, a simi- 
lar equations set depending on known values of characteristics and 
temperatures ofa previous state, k, is introduced [52]. Once temper- 
atures are known for t=k, they will be calculated for t=k+1 using 
an additional coefficient. The general coefficients from Eqs. (7.a) 
and (7.b) are changed according to Eqs. (9.a) and (9.b), respectively 
(other terms and equations are adapted accordingly). 


Tey a TEI -TEHA = aey (S+ ULTo) — rT%;; 
A= E H = Ax Ay 4 r| i ale (9.a) 
THID + TE — Ty + + Ty) 
= -ITÉ Ip = am, fori = 1andi =m (9.b) 


2.3. Solution to model optimization procedure 
The proposed procedure consists in the following steps: 


a) The system and its model are defined analytically and numer- 
ically for one and two dimensions [46,54]. Initially, the values 
of certain physical properties of the collector materials (input 
parameters of the model) are provided according to the lit- 
erature. These parameters are integrated into the vector of 
characteristics UP that may vary for each instant of time n 
(named UP’). Initially, the parameters included in this vector 
are the following: 

- ke: extinction coefficient of the covering glass (m~!) 

- Ny: refraction index of the glass 

- dg: directional absorbance of the copper black plate, consid- 
ering an incident angle of irradiation over the tilted surface, 0, 
depending on tilt, latitude, declination, solar time, etc. 

- €p: plate emissivity 


START 


PARAMETERS, UP" 
NUMERICAL 
MODELLING 


TWO-DIMENSIONAL TEMPERATURE 
DISTRIBUTION ON THE PLATE 


(Tpmodel is taken) 
ONE-DIMENSIONAL TEMPERATURE 


DISTRIBUTION OF FLOW INTO TUBES 
(Troe! is taken) 


OBJECTIVE FUNCTION, w 


MEASURED DATA: 
Tpeal , Troal 


PARAMETER 
MODIFICATION 
ALGORITHM 


Fig. 4. Flow chart of the identification process (steps a-c). 


- €g: glass emissivity 

- kir: conduction heat transfer coefficient of insulation [W/(mK)] 

- kp: conduction heat transfer coefficient of copper [W/(mK)] 

- kp: proportionality parameter associated to the tube inner wall 
film coefficient (convection heat transfer), ho [W/(m? K)] [46]. 

- kg: proportionality parameter associated to the volumetric 
equivalent heat transfer coefficient, Uc [W/(m? K)], consider- 
ing two-dimensional systems [46]. 

b) Once the system is defined, operating conditions are tested two 
hours a day. These conditions include different days and collec- 
tor tilts. For each instant of time ‘n’, the vector SP" contains the 
following input variables associated to both the irradiation and 
the system behavior (i.e. flow rate and temperatures): 

- Mm: water flow rate (kg/s) 

- s: collector tilt (rad) 

- g: latitude (rad) 

- w: timing angle (rad) 

- 6: declination (rad) 

- y: azimutal angle (rad) 

- 0: incident angle of beam irradiance over the tilted surface 
(rad) 

- Ip: global irradiance over the tilted collector (W/mz2). 

- Ta: room temperature (K) 

- Tg: fluid inlet temperature (K) 
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Fig. 5. Temperature distribution on the absorber: (a) isometric view; (b) plant view. 
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Table 2 
Characteristic angles of the collector for some random measures. 
Stationary test no. Collector tilt s (°) Latitude ø (°) Solar hour angle w (°) Declination ô (°) Azimutal angle y (°) Incident angle of direct 
irradiation 0 (°) 
1 38 38 15 20.73 0 25.56 
2 28 38 22.5 21:27 0 24.47 
3 48 38 15 21.43 0 34.9 
Table 3 
Input-output data values of random measures from Table 2. 
Stationary test Solar incident Room temperature Inlet water Outlet water Bulb plate Fluid flow rate 
no. radiation I; (W/m?) To (°C) temperature T; (°C) temperature T, (°C) temperature Tp (°C) mr (kg/h) 
1 936.8 23.2 31 50.5 57 6.42 
2 900.8 19.4 26 39 49.5 7.35 
3 862.1 25:3 33 46 50.5 7.35 
Table 4 
Optimized values for unknown parameters after the identification process shown in Tables 2 and 3. 
Test no. ke (mm7!) n2 o Ep Eg kr (W/(mK)) kp (W/(mK)) U: (W/(m? K)) h (W/(m? K)) 
Initial 0.4 1.526 0.801 0.039 0.88 0.029 402.4 11.8 x 10? 11 
Initial value for Initial value for 
test No. 1 test No. 1 
1 0.4022 1.5649 0.7543 0.0391 0.8804 0.0293 457.7 11.9 x 104 58.84 
2 0.4016 1.5667 0.7405 0.0391 0.8802 0.0292 423.71 9.1 x 104 35.7 
3 0.4019 1.5569 0.7668 0.0391 0.8803 0.0293 474.38 18.4 x 104 69 


c) To undertake the simulation, the variation of the input data of 
the vector SP” vs. time is required. On the other hand, vector 
UP", built-up with the previously mentioned unknown terms 
(provided by the literature), could lead to a weak approxima- 
tion to reality. Thus, the target of the application at this step is 
to identify the real values of those unknown terms, simulating 
the real behavior of the collector over the time. This method 
is used for steady state analysis. The objective function W is 
calculated from the sum of squares of differences between the 
model response and the empirical system response values, asso- 
ciated to both the outlet temperature of the working fluid from 
the collector Tf, (K) and the plate temperature of a represen- 
tative point of the plate, T, (K), placed near to the top of the 
absorber (the real value is measured by a bulb thermometer), as 


0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 


0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 
Transversal length (m) 


Fig. 6. Temperature variation (°C) in three cross sections of the absorber, transver- 
sally to fluid movement at different locations of the flat plate. 


shown in Figures 1 and 2. So, the objective function remains as 
follows: 

we Chau N T model f + (Teal = Tmodel y? (10) 
If hydraulic friction loss and fluid flow rate distribution are con- 
sidered, then the objective function must incorporate quadratic 
terms related to pressure in some points of the circuit. This 
additional approach has been previously evaluated in the same 
system but working as a thermosyphon, providing accurate 
results [53]. The identification procedure is applied n times to 
determine the parameters values of UP", thus minimizing the 
described function W. The identification process described in this 
step has been summarized in Fig. 4. 


d) Next step comprises the determination of SP" components (i.e. 


heat release by incident radiation and angles related to solar 


FLOW TEMPERATURE INTO RISERS 
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Fig. 7. Water temperature evolution through pipes. 
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Fig. 8. Comparison between real and predicted values for parameters kuc (a) and Nu 


(b). 


time) and the system response data (i.e. plate temperatures) 
that fit the real values of the parameters belonging to UP". If the 
previous expressions do not suffice to approximate the desired 
functions, other kind of equations must be taken into consid- 
eration, i.e. a multivariate analysis depending on the operating 
variables included in SP" [55]. Statistical fitting functions for 
these parameters must be analyzed. These functions depend 
on the operating parameters included in the SP" vector, mak- 
ing them time independent. Finally, the initial model must be 
redefined. 

e) Once the model is redefined, the transient state model is checked 
for each complete daily test. 


3. Results 


This section includes validation of the numerical model, param- 
eters identification procedure and transient regime evaluations. 
Separation between nodes in each of the 14 mentioned regions is 
calculated considering Ax= Ay=0.005 m, resulting in m=7. This 
means there are 637 nodes (7 x 91) over each absorber fin (Fig. 
2). In the plate extremes, only one riser and half absorber fin are 
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Fig. 9. Experimental and modelled results at 28° tilt: (a) water outlet temperature 
evolution; (b) absorber plate temperature evolution. 


considered. So, the number of nodes is 364 (4 x 91). Also, the fluid 
temperatures in the boundary limits, considering the x-axis, have 
been considered. Thus, the set of equations is composed of 810 
equations, with 810 unknown values of temperatures linked to each 
analyzed absorber portion. The intermediate space between nodes 
in the upper and lower tubes is 0.03 m. 


3.1. Steady state parameter identification 


Initially, the unknown vector is defined with data from the lit- 
erature. Parameters identification process of UP” is carried out 
randomly over the time. Data belonging to vectors SP" are depicted 
in Tables 2 and 3, while Table 4 shows the optimization results of 
UP" under the same conditions. 

Fig. 5 shows the temperatures distribution over the plate in a 
particular case study, while Fig. 6 exhibits three different transver- 
sal sections. In agreement with the literature, the variation is small, 
but increases with the distance from the riser [17]. The bulb ther- 
mometer shows same values for modelled and real temperatures. 
Also, extreme values of input and output water temperatures match 
real ones. Fig. 7 depicts one-dimensional temperature variation of 
the fluid flow into each riser. 

As a result of the first evaluation, it can be seen that material 
properties show values close to those initially estimated. Varia- 
tions in the glass covering, isolation and copper plate properties 
demonstrate that values from the literature are appropriate to the 
simulation process [54]. In fact, they almost satisfy the condition of 
monotonicity needed to assure an unique solution of the IHTP [40]. 
For this reason, these values are considered known and removed 
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Fig. 10. Experimental and modelled results at 38° tilt: (a) water outlet temperature 
evolution; (b) absorber plate temperature evolution. 


from the vector of unknown parameters UP". Thus, computational 
charge can be reduced significantly in the next step. On the other 
hand, parameters associated to heat transfer experiment important 
variations compared to the initial estimation (from 5 to 10 times the 
original value) provided by Eqs. (3) and (5) due to model simplifica- 
tions. The same conclusion can be drawn when analytical functions 
are applied [49]. Only two terms belonging to UP", namely kuc and 
kp, vary Uc and h values, respectively, over the time. Flow regime 
with a very low Re makes difficult to calculate the flow rate into 
each riser, the heat transfer mechanism and the real convection 
coefficient. This model assumes a two-dimension problem, consid- 
ering the pipes as flat surfaces. An accurate study of the associated 
terms Uc and h, considering physical and geometric properties, is 
required. 

Next step consists in the identification of the data collected two 
hours per day, during several days. As mentioned above, only the 
parameters kuc and kp belong to the vector UP” and will be opti- 
mized. Also, there is a relation between Uc and h, as shown in Eq. 


(5). 


3.2. New functions and model review 


Once the selected parameters are identified, a relation between 
Uc and temperature is observed. So, next step requires to find out 
a regression function for parameter kuc, providing the equivalent 
heat transfer coefficient, Uc. With the help ofa multivariate analysis, 
the parameter is adjusted by an exponential/potential function, f; 
[55]. 
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Fig. 11. Experimental and modelled results at 48° tilt: (a) water outlet temperature 
evolution; (b) absorber plate temperature evolution. 


To approximate each parameter (output variable), regression 
functions with exponential/potential shape, fı, are required. For 
each output function, fj, it can be found wl different input variables, 
Xun conforming the vector SP”. The function is as follows: 


fi = aX? XbA XO; Yue (1, w) > Xy eS (11) 


where a; and b,, are constants of the function when it is successfully 
developed. 

Initially, it is assumed that f4 depends on different temperatures, 
namely inlet and outlet fluid, room and plate temperatures. Nev- 
ertheless, this first approach is not enough to ensure a good fitting. 
So, other parameter like solar irradiation over the tilted surface is 
needed. According to this, the best result for fı [showing a standard 
deviation of 0.033 of the root mean square (RMS) error applied to 
the function] is found when the released solar heat under differ- 
ent temperatures is taken into consideration, as shown in Eq. (12) 
(where a; is constant). 


kuc = f1(@1, Tp, Th, To, To, It) (12) 


In relation with the film coefficient, h, it can be seen that final data 
present a great dispersion. This confirms the necessity of establish- 
ing a new function to approximate Nu and identify the real film 
coefficient h over the time. It shows there is a non-developed fluid 
flowing along the risers, even presenting internal recirculation, as 
indicated by Zghal [29]. So, a deep analysis of the involved pro- 
cesses, aided by CFD techniques, is required. However, this situation 
goes beyond the objectives of this article. So, a regression function 
to estimate h and kuc is enough for the proposed study. 

In the present work, kp helps to determine the real value of 
h and subsequently Nu. Then, the procedure for searching an 
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approximation function f2 by stochastic methods is repeated and 
Nuis calculated by the input variables X,,2. Between these variables, 
the non-dimensional parameters Re, Pr and Gr provide important 
information about water flowing into the risers and inner wall tem- 
peratures. Nevertheless, the simplification and the non-developed 
fluid flow could lead to misstatements during heat transmission 
modelling. So, the inclusion of other parameters is required. As a 
result, the adequate function includes the input parameters shown 
in Eq. (13) (where az is constant), besides irradiation and tilt. This 
last term affects the heat distribution and the transversal velocity 
of the fluid flow. So, determination of h into the risers can be easily 
undertaken using expression (13), as follows: 


Nuk 
h= e: 
2R; 


Results show a standard deviation of the RMS error of 0.11. Fig. 8 
depicts the correlation between the functions fı and f2. Standard 
deviations are very low, thus indicating results are in agreement 
with parameters values from Eqs. (12) and (13). 

Approaches provided by Eqs. (12) and (13) allow the validation 
of the proposed numerical model with available data. Nevertheless, 
more general and useful correlations should be considered in future 
works, thus extending the experiments to fluid flows, temperatures 
and friction losses data [24,25,28,50,51]. In this sense, increasing 
the number of sensors is crucial. Also, to improve the model, Eq. 
(3) must be substituted by Eq. (13), while kuc must be determined 
using Eq. (12) and applied to Eq. (5). 


Nu = f2(a2, Re, Pr, Gr, Tp, Thi, Tyo, To, It, O); (13) 


3.3. Transient state analyses further to model review. Final 
validations 


Once the model has been revised, the transient state is ana- 
lyzed. Considering the input data, distributions are in agreement 
with real data, as shown in Figs. 5-7, but depending on both time 
and previous state. From these distributions, two data are selected 
and compared with real values: outlet water temperature and bulb 
temperature (located on the plate). Some results are shown in Figs. 
9-11. 

Approaches present good accuracy. Although certain delays are 
observed, they are caused by the time interval. This problem can be 
fixed by increasing the computing time, but it does not guarantee 
better approximations. In any case, the previous model, in addi- 
tion to new approximation functions ensures enough accuracy to 
simulate the real behavior of this particular system. 

Similar results are found by other researchers, considering the 
transient behavior of the collector [35], as shown in Figs. 9-11. 
Although results validate the proposed model, a deeper study of the 
thermo-hydraulic mechanisms may be considered in future works. 
Nevertheless, it is important to notice that the target of this study 
was to estimate the temperature distribution over the plate and 
into de pipes, and this objective has been successfully reached. 


4. Conclusions 


A general strategy that helps to validate a solar collector model 
where some parameters are difficult to determine, under both non- 
steady and non-uniform states, has been proposed. In this case, the 
specific conditions of the fluid regime and the geometry of the col- 
lector make difficult to understand the dynamic behavior of the 
system. To avoid the use of complex numerical techniques, such 
as CFD or finite element method (FEM), a highly accurate model 
including simplifications of the physical unit (plate-pipes) has been 
developed. As an advantage, this process only requires a few sen- 
sors. For this reason, it makes possible the analyses in equipments 
without many sophisticated sensors (as used in heat-pipes stud- 
ies). Temperatures distribution of both absorber and water gives 


qualitatively similar results to previous works found in the litera- 
ture. Moreover, both steady and transient regimes can be validated 
with experimental results, thus demonstrating the accuracy of the 
model. Results also show that materials properties are in close 
agreement with the values given by the literature. This means low 
uncertainties if similar values are used, thus avoiding the need of 
highly accurate initial data. 

Apart from the type of model, both the identification process and 
the approximation functions constitute two fundamental steps. 
Nevertheless, the processes depend on the following conditions: 


1. The identification process is based on a non-linear optimization 
procedure over a quadratic continuous function that is a deriva- 
tive dependent function. In case the objective function exhibits 
a high discontinuity, the proposed methodology could be inef- 
fective. If so, to optimize the function, other techniques such as 
Genetic Algorithms must be evaluated. 

2. Considering parameters approximation functions, exponential 
regression functions have been developed by multivariate anal- 
ysis of variance. However, other functions may also be evaluated, 
i.e. linear functions and potential functions. Including some 
modifications to the exponential regression functions provides 
another possibility to take into consideration. In this case, 
stochastic models cannot be used. So, the aid of CFD techniques 
is recommended. 


The proposed methodology allows the analysis of the system under 
a huge range of different conditions, such as materials and geome- 
try changes. Also, it makes possible to find out the optimum value of 
the pipes diameter and the distance between them to maximize the 
captured energy, among many other possibilities. A further analy- 
sis to optimize the geometry and dimensions may be developed. In 
this sense, a future target could be to modify them to maximize the 
collection of energy, as mentioned in the introduction of this paper. 
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